subroutine meshCellOut_def_old(xg,n,meshCellOut,nCellsOut,i,tol)
implicit none
	real *8 xg(3)
	integer nCellsOut,i
	real *8 n,tol
	integer meshCellOut(*)
	real *8 d,rr1,rr2,vtmp
	vtmp=0.14d0*n
	rr1=vtmp-tol
	rr2=vtmp+tol
	if (xg(1).gt.rr1.and.xg(1).lt.rr2 ) then
	    if (dabs(xg(2)).lt.10d-4.or.dabs(xg(3)).lt.10d-4) then
            nCellsOut=nCellsOut+1
            meshCellOut(nCellsOut)=i
    	    return
        end if
	    if (dabs(xg(3)).lt.10d-4.or.dabs(xg(2)).lt.10d-4) then
            nCellsOut=nCellsOut+1
            meshCellOut(nCellsOut)=i
        end if
    end if
	
	return
end subroutine meshCellOut_def_old

subroutine meshCellOut_def(xg,n,meshCellOut,nCellsOut,i,tol)
implicit none
	real *8 xg(3)
	integer nCellsOut,i
	real *8 n,tol,dtol
	integer meshCellOut(*)
	real *8 d,rr1,rr2,vtmp

	vtmp=0.14d0*n
	dtol=5d-3
	rr1=vtmp-tol
	rr2=vtmp+tol
	if (xg(1).gt.rr1.and.xg(1).lt.rr2 ) then
	    if (dabs(xg(2)).lt.dtol.or.dabs(xg(3)).lt.dtol) then
            nCellsOut=nCellsOut+1
            meshCellOut(nCellsOut)=i
	        return
        end if
	    if (dabs(xg(3)).lt.dtol.or.dabs(xg(2)).lt.dtol) then
            nCellsOut=nCellsOut+1
            meshCellOut(nCellsOut)=i
        end if
    end if
	return
end subroutine meshCellOut_def

subroutine meshCellOut_defn(x,meshCellOut,nCellsOut,iel, &
    	area,xCenter,meshSf, nUnusedSections)
!           meshCellOut_defn(xt,meshCellOut,nCellsOut,i,&
!		areaOut(nCellsOut+1),xOut(1:3,nCellsOut+1),normalOut(1:3,nCellsOut+1))
!use Data_cfoam, only: AcousticsOnOff
use Data_cfoam, only: ZeroEpsilon
use Acoustics, only: AcousticSurf, DefNormalDirection
implicit none
	real(8),intent(in) :: x(3,8)
	integer,intent(in) :: iel
	real(8),intent(out) :: area,xCenter(3),meshSf(3)
	integer,intent(inout) :: nCellsOut, nUnusedSections
	integer,intent(inout) :: meshCellOut(*)
	
	integer rib(2,12)
	real *8 xrib(3,12),xCells(3,12)
	integer location,i,j,k,locationp,nPoints,nVar(4,6),n1,n2,n3,n4
	real *8 intersect(12),mag,z(3), mag1
!	real(8) :: r(3)
	
	
	data (rib(1,i),i=1,12)/1,2,3,4,1,2,3,4,5,6,7,8/
	data (rib(2,i),i=1,12)/2,3,4,1,5,6,7,8,6,7,8,5/
	data (nVar(1,i),i=1,6)/1,1,1,1,1,1/
	data (nVar(2,i),i=1,6)/2,2,3,3,4,4/
	data (nVar(3,i),i=1,6)/3,4,2,4,2,3/
	data (nVar(4,i),i=1,6)/4,3,4,2,3,2/
	
!	write(*,*) x
	call surface(x(1:3,1),location)
	if(location.eq.0) return
	do i=2,8
	    call surface(x(1:3,i),locationp)
	    if(locationp.eq.0) return
	end do
	do i=2,8
	    locationp=location
	    call surface(x(1:3,i),location)
!	    if(locationp.eq.0) return
	    if (locationp.ne.location) then
!	write(*,*) nCellsOut+1,iel
            nCellsOut=nCellsOut+1
            meshCellOut(nCellsOut)=iel
	        nPoints=0
	        xCenter=0d0
	        do j=1,12
	            call line(x(1:3,rib(1,j)),x(1:3,rib(2,j)),xrib(1:3,j),intersect(j))
	            if (intersect(j).ge.0d0) then
	                nPoints=nPoints+1
	                xCells(1:3,nPoints)=xrib(1:3,j)
	                xCenter(1:3)=xCenter(1:3)+xrib(1:3,j)
!	                write(*,*) j,intersect(j)
	            end if
	        end do
	        xCenter=xCenter/nPoints

!	        if (nPoints.gt.4) stop 'Fatal Error in FW'
            if (nPoints.gt.4)then
                write(*,*) '@Warning: nPoints.gt.4'
                nCellsOut=nCellsOut-1
                nUnusedSections = nUnusedSections + 1
            end if
	        if (nPoints.eq.3) then

                call normal2TriangleDef(meshSf(1:3),xCells(1:3,1),xCells(1:3,2),xCells(1:3,3))
                
                area = sqrt(sum(meshSf(1:3)*meshSf(1:3)))
                if(area.lt.ZeroEpsilon)then
                    meshSf(1) = 1.0
                    meshSf(2:3) = 0.0
                    area = 0d0
                    return
                end if

!                call sxdef(meshSf(1:3),xCells(1:3,1),xCells(1:3,2),xCells(1:3,3),xCells(1:3,3)) !,xCenter(1:3),xCenter(1:3))
!	            mag=dsqrt(meshSf(1)**2+meshSf(2)**2+meshSf(3)**2)
!	            if(mag.lt.1d-8)then
!	                meshSf(1) = 1.0
!	                meshSf(2:3) = 0.0
!	                area = 0.0
!	                return
!	            end if

!	            meshSf=meshSf/mag
                
                meshSf = meshSf / area
                if(AcousticSurf)then ! AcousticsOnOff
                    call DefNormalDirection(xCenter,meshSf)
                else
	                z=xCenter+meshSf *1d-1
	                call surface(z,location)
	                if (location.lt.0) meshSf=-meshSf
                end if

!                call sxdef(z,xCells(1:3,1),xCells(1:3,2),xCells(1:3,3),xCells(1:3,3))
!	            area=dsqrt(z(1)**2+z(2)**2+z(3)**2)
	            
	            return
	        end if

! ***************************************************************************************** !
	        call CellOrderDef(xCells(1:3,1:4))

	        xCenter(1:3) = 0d0
	        do j = 1,4
    	        xCenter(1:3) = xCenter(1:3) + xCells(1:3,j)
    	    end do
            xCenter(1:3) = 0.25d0 * xCenter(1:3)

            meshSf = 0d0
            do j = 1,4
!                call sxdef(z(1:3),xCells(1:3,j),xCells(1:3,j+1),xCenter(1:3),xCenter(1:3))
                call normal2TriangleDef(z(1:3),xCells(1:3,j),xCells(1:3,mod(j,4)+1),xCenter(1:3))
                meshSf(1:3) = meshSf(1:3) + z(1:3)
            end do
            area = sqrt(sum(meshSf(1:3)*meshSf(1:3)))
            if(area.lt.ZeroEpsilon)then
                meshSf(1) = 1.0
                meshSf(2:3) = 0.0
                area = 0d0
                return
            end if

            meshSf = meshSf / area
	        if(AcousticSurf)then ! AcousticsOnOff)then
                call DefNormalDirection(xCenter,meshSf)
            else
	            z=xCenter+meshSf *1d-1
	            call surface(z,location)
	            if (location.lt.0) meshSf=-meshSf
            end if

! ***************************************************************************************** !
!	        call sxdef(meshSf(1:3),xCells(1:3,1),xCells(1:3,2),xCells(1:3,3),xCells(1:3,3))
!	        mag=dsqrt(meshSf(1)**2+meshSf(2)**2+meshSf(3)**2)
!	        
!	        call sxdef(z(1:3),xCells(1:3,1),xCells(1:3,4),xCells(1:3,3),xCells(1:3,3))
!	        mag1=dsqrt(z(1)**2+z(2)**2+z(3)**2)
!	        
!	        if(mag.lt.mag1)then
!	            meshSf(1:3) = z(1:3)
!	            mag = mag1
!	        end if
!	        
!	        if(mag.lt.1d-8)then
!                meshSf(1) = 1.0
!                meshSf(2:3) = 0.0
!                area = 0.0
!                return
!            end if
!	        meshSf=meshSf/mag
!	        if(AcousticSurf)then ! AcousticsOnOff)then
!                call DefNormalDirection(xCenter,meshSf)
!            else
!	            z=xCenter+meshSf *1d-1
!	            call surface(z,location)
!	            if (location.lt.0) meshSf=-meshSf
!            end if
!	        
!	        area=0d0
!	        do k=1,6
!	            n1=nVar(1,k)
!	            n2=nVar(2,k)
!	            n3=nVar(3,k)
!	            n4=nVar(4,k)
!                call sxdef(z,xCells(1:3,n1),xCells(1:3,n2),xCells(1:3,n3),xCells(1:3,n4))
!	            mag=dsqrt(z(1)**2+z(2)**2+z(3)**2)
!	            if (area.lt.mag) area=mag
!	        end do
! ***************************************************************************************** !

!    	    r(1:3) = xCenter(1:3) + meshSf(1:3)


    !	stop
	        return
	    end if
	end do

	return
end subroutine meshCellOut_defn
	
subroutine surface(x,n)
use Acoustics, only: AcousticSurf, HeavysideFunc
implicit none
	real(8),intent(in) :: x(3)
	integer,intent(out) :: n
	
	real *8 rr,rr1
	
	if(AcousticSurf)then
	    n = HeavysideFunc(x)
	    return
	end if
	
	n=1
	rr=dsqrt(x(2)**2+x(3)**2)
	rr1=0.17 !*28
!	rr1=0.158+0.158*14.0/180.0*3.1415927*x(1)
	if (rr.lt.rr1) n=-1
	return
end subroutine surface

subroutine line(x,y,z,alpha)
implicit none
	real(8),intent(in) :: x(3),y(3)
	real(8),intent(out) :: z(3),alpha
	
	real *8 x0,x1,x2
	integer loc,loc1,loc2
	
	x1=0d0
	x2=1d0
	call surface(x,loc1)
	call surface(y,loc2)
	alpha=-1d0
	if (loc1.eq.loc2) return
 1	x0=(x1+x2)/2d0
 	z(1:3)=x(1:3)*(1.0d0-x0)+y(1:3)*x0
	call surface(z,loc)
	if (loc.eq.loc1) x1=x0
	if (loc.eq.loc2) x2=x0
	if (dabs(x1-x2).gt.1d-6) goto 1
	alpha=x0
	
	return
end subroutine line


subroutine meshCellTransducerOut(x,iel,Hollow,insideOnly,nCellsOut,meshCellOut)
use Acoustics, only: HeavysideFunc
implicit none
	real(8),intent(in) :: x(3,8)
	integer,intent(in) :: iel
	logical,intent(in) :: Hollow, insideOnly
	integer,intent(inout) :: nCellsOut
	integer,intent(inout) :: meshCellOut(*)
	
	integer rib(2,12)
	integer location,i,locationp,nVar(4,6)

	data (rib(1,i),i=1,12)/1,2,3,4,1,2,3,4,5,6,7,8/
	data (rib(2,i),i=1,12)/2,3,4,1,5,6,7,8,6,7,8,5/
	data (nVar(1,i),i=1,6)/1,1,1,1,1,1/
	data (nVar(2,i),i=1,6)/2,2,3,3,4,4/
	data (nVar(3,i),i=1,6)/3,4,2,4,2,3/
	data (nVar(4,i),i=1,6)/4,3,4,2,3,2/
	
	if(Hollow)then
	    location = HeavysideFunc(x(1:3,1))
	    if(location.eq.0) return
	    do i=2,8
	        locationp = HeavysideFunc(x(1:3,i))
	        if(locationp.eq.0) return
	    end do
	    do i=2,8
	        locationp = HeavysideFunc(x(1:3,i))

	        if (locationp.ne.location) then
                nCellsOut=nCellsOut+1
                meshCellOut(nCellsOut)=iel
                return
            end if
        end do
    else
	    do i=1,8
	        location = HeavysideFunc(x(1:3,i))

            if(insideOnly)then
                if(location.ne.-1) return
            else
    	        if (location.eq.-1) then
                    nCellsOut=nCellsOut+1
                    meshCellOut(nCellsOut)=iel
                    return
                end if
            end if
        end do
        
        if(insideOnly)then
            nCellsOut=nCellsOut+1
            meshCellOut(nCellsOut)=iel
        end if

    end if

end subroutine meshCellTransducerOut

subroutine CellOrderDef(xCells)
implicit none
    real(8),intent(inout) :: xCells(1:3,1:4)
    
    real(8) :: L(3,3), a,b, r(3)
    integer :: i,j
    real(8) :: Laa,Lab,Lbb, Lac,Lbc, D
    
    do j = 1,3
        do i = 1,3
            L(i,j) = xCells(i,j+1) - xCells(i,1)
        end do
    end do

    Laa = sum(L(1:3,1) * L(1:3,1))
    Lab = sum(L(1:3,1) * L(1:3,2))
    Lbb = sum(L(1:3,2) * L(1:3,2))

    D = Laa*Lbb - Lab**2

    Lac = sum(L(1:3,1) * L(1:3,3))
    Lbc = sum(L(1:3,2) * L(1:3,3))

    a = ( Lbb * Lac - Lab * Lbc) / D
    b = (-Lab * Lac + Laa * Lbc) / D

    if(a.gt.0.and.b.gt.0)then
        call ExchangeVectors(xCells(1:3,3),xCells(1:3,4))
        return
    end if

    if(a.lt.0.and.b.gt.0)then
        return
    end if

    if(a.lt.0.and.b.lt.0)then
        ! degenerate case
        xCells = 0d0
        return
    end if

    if(a.gt.0.and.b.lt.0)then
        call ExchangeVectors(xCells(1:3,1),xCells(1:3,4))
        return
    end if

end subroutine CellOrderDef

subroutine ExchangeVectors(a,b)
implicit none
    real(8),intent(inout) :: a(3),b(3)
    real(8) :: r(3)
    
    r(1:3) = a(1:3)
    a(1:3) = b(1:3)
    b(1:3) = r(1:3)

end subroutine ExchangeVectors

subroutine normal2TriangleDef(z,a,b,c)
implicit none
    real(8),intent(in) :: a(3),b(3),c(3)
    real(8),intent(out) :: z(3)

    real(8) :: l1(3),l2(3)


    l1 = a - c
    l2 = b - c
    call VectorProduct(l1,l2,z)

    z =  0.5d0 * z

end subroutine normal2TriangleDef

subroutine VectorProduct(l1,l2,l)
implicit none
    real(8),intent(in) :: l1(3),l2(3)
    real(8),intent(out) :: l(3)

    l(1) = l1(2)*l2(3) - l1(3)*l2(2)
    l(2) = l1(3)*l2(1) - l1(1)*l2(3)
    l(3) = l1(1)*l2(2) - l1(2)*l2(1)

end subroutine VectorProduct